* MASTER_FIGURE.do
clear all
graph set window fontface "Ariel"

* *******************************************************
* FIGURE F.2 IN APPENDIX -- IDENTIFICATION FIGURES ******
* *******************************************************
import delimited "../../output/quant/output_calib_identidicationtest_change.csv"

#d;
twoway (scatter db_base cbddist if cbddist<=3, sort color(black) msize(small) )
(lpoly db_base cbddist if cbddist<=3, sort lcolor(black) degree(6) lwidth(thick))
(scatter db_weak cbddist if cbddist<=3, sort color(navy) msymbol(Oh) msize(small))
(lpoly db_weak cbddist if cbddist<=3, sort  lcolor(navy) lpattern(dash) degree(6) lwidth(thick))
(scatter db_strong cbddist if cbddist<=3, sort color(cranberry) msymbol(Sh) msize(small))
(lpoly db_strong cbddist if cbddist<=3, sort lpattern(shortdash) lcolor(cranberry) degree(6) lwidth(thick)), 
	legend(order(1 "Estimated agglomeration" 3 "Weak agglomeration" 5 "Strong agglomeration" 2 "Estimated agglomeration (fitted line)" 4 "Weak agglomeration (fitted line)" 6 "Strong agglomeration (fitted line)") row(2) size(small))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.55)
	xtitle("Distance to CBD (kilometers)", size(5.5)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Average change in log value of idiosyncratic term" "of fundamental amenity in 1955-75", size(5.5)) ylabel( -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08,angle(horizontal) labsize(large)) ;
	graph export "../../output/figure/figure_identification_amenity_change.pdf", as(pdf)  replace ;
#d cr

#d;
twoway  (scatter da_base cbddist if cbddist<=3, sort color(black) msize(small))
(lpoly da_base cbddist if cbddist<=3, sort lcolor(black) degree(6) lwidth(thick))
(scatter da_weak cbddist if cbddist<=3, sort color(navy) msymbol(Oh) msize(small))
(lpoly da_weak cbddist if cbddist<=3, sort  lpattern(dash) lcolor(navy) degree(6) lwidth(thick))
(scatter da_strong cbddist if cbddist<=3, sort color(cranberry) msymbol(Sh) msize(small))
(lpoly da_strong cbddist if cbddist<=3, sort  lpattern(dash) lcolor(cranberry) degree(6) lwidth(thick)), 
	legend(order(1 "Estimated agglomeration" 3 "Weak agglomeration" 5 "Strong agglomeration" 2 "Estimated agglomeration (fitted line)" 4 "Weak agglomeration (fitted line)" 6 "Strong agglomeration (fitted line)") row(2) size(small))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.55)
	xtitle("Distance to CBD (kilometers)", size(5.5)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Average change in log value of idiosyncratic term" "of fundamental productivity in 1955-75", size(5.5)) ylabel( -0.06 -0.04 -0.02 0 0.02 0.04 0.06,angle(horizontal) labsize(large)) ;
	graph export "../../output/figure/figure_identification_productivity_change.pdf", as(pdf)  replace ;
#d cr

* ***********************************************************
* FIGURE F.3 IN APPENDIX -- IDENTIFICATION FIGURES **********
* ***********************************************************
import delimited "../../output/quant/output_calib_identidicationtest_level.csv", clear
gen l_area = ln(area)
reg b_base l_area 
predict b_base_narea, residual 
reg b_weak l_area 
predict b_weak_narea, residual 
reg b_strong l_area 
predict b_strong_narea, residual 

reg a_base l_area 
predict a_base_narea, residual 
reg a_weak l_area 
predict a_weak_narea, residual 
reg a_strong l_area 
predict a_strong_narea, residual 

#d;
twoway (lpoly b_base_narea cbddist if cbddist<=3, sort lcolor(black) degree(6) lwidth(thick))
(lpoly b_weak_narea cbddist if cbddist<=3, sort  lcolor(navy) lpattern(dash) degree(6) lwidth(thick))
(lpoly b_strong_narea cbddist if cbddist<=3, sort lpattern(shortdash) lcolor(cranberry) degree(6) lwidth(thick)), 
	legend(order(1 "Estimated agglomeration" 2 "Weak agglomeration" 3 "Strong agglomeration") row(1) size(medium))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.55)
	xtitle("Distance to CBD (kilometers)", size(5.5)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Average log value of idiosyncratic term" "of fundamental amenity in 1955-75 (net area size of block)", size(5.5)) ylabel(-0.1 -0.05 0 0.05 0.1,angle(horizontal) labsize(large)) ;
	graph export "../../output/figure/figure_identification_amenity_level.pdf", as(pdf)  replace ;
#d cr

#d;
twoway (lpoly a_base_narea cbddist if cbddist<=3, sort lcolor(black) degree(6) lwidth(thick))
(lpoly a_weak_narea cbddist if cbddist<=3, sort  lcolor(navy) lpattern(dash) degree(6) lwidth(thick))
(lpoly a_strong_narea cbddist if cbddist<=3, sort lpattern(shortdash) lcolor(cranberry) degree(6) lwidth(thick)), 
	legend(order(1 "Estimated agglomeration" 2 "Weak agglomeration" 3 "Strong agglomeration") row(1) size(medium))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.55)
	xtitle("Distance to CBD (kilometers)", size(5.5)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Average log value of idiosyncratic term" "of fundamental productivity in 1955-75 (net area size of block)", size(5.5)) ylabel(-0.2 -0.15 -0.1 -0.05 0 0.05 0.1,angle(horizontal) labsize(large)) ;
	graph export "../../output/figure/figure_identification_productivity_level.pdf", as(pdf)  replace ;
#d cr

* CALIBRATED AMENITIES AND PRODUCTIVITY 
import delimited "../../output/quant/output_calib_fundamentals.csv", clear
save temp, replace 

import delimited "../../output/quant/output_fundamentals_1945_50.csv", clear
save temp_main, replace 

merge 1:1 geocode using "temp.dta"
drop _merge 

// Log Change in Idiosyncratic Terms from 1955 to 75 
gen bbdif60=ln(bb_err60/bb_err55)
gen aadif60=ln(aa_err60/aa_err55)
 
gen bbdif65=ln(bb_err65/bb_err60)
gen aadif65=ln(aa_err65/aa_err60) 

gen bbdif70=ln(bb_err70/bb_err65)
gen aadif70=ln(aa_err70/aa_err65) 

gen bbdif75=ln(bb_err75/bb_err70)
gen aadif75=ln(aa_err75/aa_err70)

// Log Change in Idiosyncratic Terms from 1950 to 55 
gen bbdif55=ln(bb_err55/bb_err50)
gen aadif55=ln(aa_err55/aa_err50)

// normalize fundamental amenity(b) and productivity(a) 
foreach i in 55 60 65 70 75{
	gen l_bb`i' = ln(bb`i')
	egen l_bb`i'_ave = mean(l_bb`i')
	gen l_bb`i'_v = l_bb`i' - l_bb`i'_ave 
	gen bb`i'_v = exp(l_bb`i'_v) 
	gen l_aa`i' = ln(aa`i')
	egen l_aa`i'_ave = mean(l_aa`i')
	gen l_aa`i'_v = l_aa`i' - l_aa`i'_ave 
	gen aa`i'_v = exp(l_aa`i'_v) 
}

gen l_aa_ave_55_75 = (l_aa55 + l_aa60 + l_aa65 + l_aa70 + l_aa75)/5 
gen l_bb_ave_55_75 = (l_bb55 + l_bb60 + l_bb65 + l_bb70 + l_bb75)/5 

// fundamentals for 1950 
gen l_bb50 = ln(bb50)
gen l_aa50 = ln(aa50)
// log area
gen l_area = log(area) 

// Fundamental Amenity and Productivity every 5 years 
gen l_aa_ave = (l_aa50 + l_aa55 + l_aa60 + l_aa65 + l_aa70 + l_aa75)/6
gen l_bb_ave = (l_bb50 + l_bb55 + l_bb60 + l_bb65 + l_bb70 + l_bb75)/6

// CBD distance and fundamental productivity (Log, area size adjusted by regression) 
foreach i in 50 55 60 65 70 75{
	reg l_bb`i' l_area 
	predict l_bb`i'_narea, residual 
	
	reg l_aa`i' l_area 
	predict l_aa`i'_narea, residual
}

reg l_aa_ave_55_75 l_area 
predict l_aa_ave_narea, residual 
reg l_bb_ave_55_75 l_area 
predict l_bb_ave_narea, residual 

* ***********************************************
* FIGURE F.7 IN APPENDIX ************************
* ***********************************************
* CBD distance and fundamental amenity (FIGURE F.7 IN APPENDIX)
#d;
twoway  (scatter l_bb50_narea cbddist, sort color(cranberry) msize(small))
(scatter l_bb55_narea cbddist, sort msymbol(Sh) color(gray) msize(small))
(scatter l_bb60_narea cbddist, sort msymbol(Sh) color(black) msize(small))
(scatter l_bb65_narea cbddist, sort msymbol(Oh) color(gray) msize(small))
(scatter l_bb70_narea cbddist, sort msymbol(Oh) color(black) msize(small))
(scatter l_bb75_narea cbddist, sort msymbol(Th) color(gray) msize(small))
(lpolyci l_bb50_narea cbddist, level(95) degree(6) nofit fcolor(gray%20)), 
	legend(order(1 "1950" 2 "1955" 3 "1960" 4 "1965" 5 "1970" 6 "1975" 7 "95% CI for 1950") row(1) size(medsmall))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (kilometers)", size(5.5)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Log value of fundamental amenities (b)" "adjusted with area size", size(5.5)) ylabel(-0.4 -0.2 0 0.2 0.4,angle(horizontal) labsize(large)) yscale(range(-0.4 0.4)) ;
	graph export "../../output/figure/figure_fundamental_amenity_allperiods.pdf", as(pdf)  replace ;
#d cr

* CBD distance and fundamental productivity (FIGURE F.7 IN APPENDIX)
#d;
twoway  (scatter l_aa50_narea cbddist, sort color(cranberry) msize(small))
(scatter l_aa55_narea cbddist, sort msymbol(Sh) color(gray) msize(small))
(scatter l_aa60_narea cbddist, sort msymbol(Sh) color(black) msize(small))
(scatter l_aa65_narea cbddist, sort msymbol(Oh) color(gray) msize(small))
(scatter l_aa70_narea cbddist, sort msymbol(Oh) color(black) msize(small))
(scatter l_aa75_narea cbddist, sort msymbol(Th) color(gray) msize(small))
(lpolyci l_aa50_narea cbddist, level(95) degree(6) nofit fcolor(gray%20)),
	legend(order(1 "1950" 2 "1955" 3 "1960" 4 "1965" 5 "1970" 6 "1975" 7 "95% CI for 1950") row(1) size(medsmall))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (kilometers)", size(5.5)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Log value of fundamental productivity (a)" "adjusted with area size", size(5.5)) ylabel(-0.4 -0.2 0 0.2 0.4,angle(horizontal) labsize(large)) yscale(range(-0.4 0.4)) ;
	graph export "../../output/figure/figure_fundamental_productivity_allperiods.pdf", as(pdf)  replace ;
#d cr

* *************************************************
* FIGURE F.4 IN APPENDIX  *************************
* ************************************************* 
#d;
twoway (lpoly l_aa55_narea cbddist [aw=area], sort degree(3) lwidth(thick) lpattern(dash) color(gray))
(lpoly l_aa65_narea cbddist [aw=area], sort degree(3) lwidth(thick) lpattern(longdash) color(black))
(lpoly l_aa75_narea cbddist [aw=area], sort degree(3) lwidth(thick) lpattern(line) color(cranberry))
(scatter l_aa_ave_narea cbddist, sort color(black)), 
	legend(order(1 "1955" 2 "1965" 3 "1975" 4 "Average between 1955-75") row(1) size(large))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (kilometers)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log value of fundamental productivity" "adjusted with size", size(7.0)) ylabel(-0.2 -0.1 0 0.1 0.2, angle(horizontal) labsize(huge)) yscale(range(-0.25 0.25)) ;
	graph export "../../output/figure/figure_fundamental_log_productivity_netarea.pdf", as(pdf)  replace ;
#d cr

#d;
twoway (lpoly l_bb55_narea cbddist [aw=area], sort degree(3) lwidth(thick) lpattern(dash) color(gray))
(lpoly l_bb65_narea cbddist [aw=area], sort degree(3) lwidth(thick) lpattern(longdash) color(black))
(lpoly l_bb75_narea cbddist [aw=area], sort degree(3) lwidth(thick) lpattern(line) color(cranberry))
(scatter l_bb_ave_narea cbddist, sort color(black)), 
	legend(order(1 "1955" 2 "1965" 3 "1975" 4 "Average between 1955-75") row(1) size(large))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (kilometers)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log value of fundamental amenities" "adjusted with area size", size(7.0)) ylabel(-0.2 -0.1 0 0.1 0.2, angle(horizontal) labsize(huge)) yscale(range(-0.25 0.25)) ;
	graph export "../../output/figure/figure_fundamental_log_amenity_netarea.pdf", as(pdf)  replace ;
#d cr

********************************************************************************
* FIGURE F.5 IN APPENDIX: CBD distance and residual amenities change *********** 
* ******************************************************************************

#d;
twoway (scatter bbdif60 cbddist, sort color(black%30))
(lpoly bbdif60 cbddist, sort degree(3) lwidth(thick) lpattern(line) color(black))
(scatter bbdif65 cbddist, sort color(gray) msymbol(triangle) color(black%50))
(lpoly bbdif65 cbddist, sort degree(3) lwidth(thick) lpattern(longdash) color(black%80))
(scatter bbdif70 cbddist , sort color(black) msymbol(square) color(cranberry%50))
(lpoly bbdif70 cbddist, sort degree(3) lwidth(thick) lpattern(dash) color(cranberry))
(scatter bbdif75 cbddist , sort color(black) msymbol(+) color(red%50))
(lpoly bbdif75 cbddist, sort degree(3) lwidth(thick) lpattern(shortdash) color(red)),
	legend(order(1 "1955-60" 2  "1955-60 (fitted line)" 3 "1960-65" 4 "1960-65 (fitted line)" 5 "1965-70" 6 "1970-75 (fitted line)" 7 "1970-75" 8 "1970-75 (fitted line)") row(2) size(large))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.5)
	xtitle("Distance to CBD (kilometers)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Change in log of residual amenities", size(7.0)) ylabel(-0.4 -0.2 0 0.2 0.4, angle(horizontal) labsize(huge))  yscale(range(-0.4 0.4)) ;
	graph export "../../output/figure/figure_fundamental_amenity_momentcondition.pdf", as(pdf)  replace ;
#d cr

********************************************************************************
* FIGURE F.6 IN APPENDIX: CBD distance and residual productivity change ********
* ******************************************************************************

#d;
twoway (scatter aadif60 cbddist, sort color(black%30))
(lpoly aadif60 cbddist, sort degree(3) lwidth(thick) lpattern(line) color(black))
(scatter aadif65 cbddist, sort color(gray) msymbol(triangle) color(black%50))
(lpoly aadif65 cbddist, sort degree(3) lwidth(thick) lpattern(longdash) color(black%80))
(scatter aadif70 cbddist , sort color(black) msymbol(square) color(cranberry%50))
(lpoly aadif70 cbddist, sort degree(3) lwidth(thick) lpattern(dash) color(cranberry))
(scatter aadif75 cbddist , sort color(black) msymbol(+) color(red%50))
(lpoly aadif75 cbddist, sort degree(3) lwidth(thick) lpattern(shortdash) color(red)),
	legend(order(1 "1955-60" 2  "1955-60 (fitted line)" 3 "1960-65" 4 "1960-65 (fitted line)" 5 "1965-70" 6 "1970-75 (fitted line)" 7 "1970-75" 8 "1970-75 (fitted line)") row(2) size(large))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.5)
	xtitle("Distance to CBD (kilometers)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Change in log of residual amenities", size(7.0)) ylabel(-0.4 -0.2 0 0.2 0.4, angle(horizontal) labsize(huge))  yscale(range(-0.4 0.4)) ;
	graph export "../../output/figure/figure_fundamental_productivity_momentcondition.pdf", as(pdf)  replace ;
#d cr

save temp.dta, replace 

* import fundamentals in 1950 for evaluation 
import delimited "../../output/quant/output_evaluations_1945_50.csv", clear 
merge 1:1 geocode using "temp.dta"

gen lpopdens45 = log(pop45/area)
gen lpopdens50 = log(pop50/area)
gen lpopdens50_p = log(pop50_pred/area)
gen lempdens45 = log(emp45/area) 
gen lempdens50 = log(emp50/area)
gen lempdens50_p = log(emp50_pred/area)

* ************************************************************************
* FIGURE 5 IN PAPER: population/employment in 1950 for evaluation ******** 
* ************************************************************************
* Log population density -- Data vs Simulation ** Add weight by area size for polyfit 
#d;
twoway (lpoly lpopdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lpopdens50_p cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lpopdens50 cbddist [aw=area], sort color(gray)lwidth(thick) lpattern(dash))
(scatter lpopdens50_p cbddist, sort color(black)) 
(scatter lpopdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log population density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) yscale(range(2 12)); 
	graph export "../../output/figure/figure_density_pop_baseline_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr

* Log Employment density -- Data vs Simulation ** Add weight by area size for polyfit
#d;
twoway (lpoly lempdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lempdens50_p cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lempdens50 cbddist [aw=area], sort color(gray) lwidth(thick) lpattern(dash))
(scatter lempdens50_p cbddist , sort color(black) ) 
(scatter lempdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log employment density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) ; 
	graph export "../../output/figure/figure_density_emp_baseline_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr


* ***********************************************************************
* FIGURE F.8 IN APPENDIX -- import 1930s results ************************
* ***********************************************************************
import delimited "../../output/quant/output_calib_fundamentals_1930s.csv", clear 

gen l_area = ln(area)
gen l_aa30 = ln(aa30)
gen l_bb30 = ln(bb30)
gen l_aa55 = ln(aa55)
gen l_bb55 = ln(bb55)

* adjust by area size 
reg l_aa30 l_area 
predict l_aa30_narea, residual

reg l_bb30 l_area 
predict l_bb30_narea, residual

reg l_aa55 l_area 
predict l_aa55_narea, residual

reg l_bb55 l_area 
predict l_bb55_narea, residual

#d;
twoway(scatter l_bb55_narea cbddist , sort color(gray) msymbol(square) msize(small))
(lpoly l_bb55_narea cbddist, lwidth(medthick) color(cranberry))
(scatter l_bb30_narea cbddist, sort color(black)  msize(small))
(lpoly l_bb30_narea cbddist, lwidth(medthick) lpattern(dash) color(cranberry)), 
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.8)
	legend(order(1 "1955" 2 "1955 fitted line" 3 "1936" 4 "1936 fitted line") row(1))
	xtitle("Distance to CBD (kilometers)", size(5)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Log value of fundamental amenities" "adjusted with area size", size(5)) ylabel(-0.2 -0.1 0 0.1 0.2, angle(horizontal) labsize(large)) yscale(range(-0.25 0.25));
	graph export "../../output/figure/figure_fundamental_amenity_1930s_cbddist.pdf", as(pdf)  replace ;
#d cr

#d;
twoway(scatter l_aa55_narea cbddist , sort color(gray) msymbol(square) msize(small))
(lpoly l_aa55_narea cbddist, lwidth(medthick) color(cranberry))
(scatter l_aa30_narea cbddist, sort color(black)  msize(small))
(lpoly l_aa30_narea cbddist, lwidth(medthick) lpattern(dash) color(cranberry)), 
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.8)
	legend(order(1 "1955" 2 "1955 fitted line" 3 "1938" 4 "1938 fitted line") row(1))
	xtitle("Distance to CBD (kilometers)", size(5)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Log value of fundamental productivity" "adjusted with area size", size(5)) ylabel(-0.2 -0.1 0 0.1 0.2,angle(horizontal) labsize(large)) yscale(range(-0.25 0.25));
	graph export "../../output/figure/figure_fundamental_productivity_1930s_cbddist.pdf", as(pdf)  replace ;
#d cr

* ***********************************************************
* FIGURE G.4 IN APPENDIX -- FUNDAMENTAL_1950_NOAGG_FIGURE ***
* ***********************************************************
import delimited "../../output/quant/output_noaggfund_50.csv", clear 

gen l_area = ln(area)
gen l_a50noagg = ln(y50noagg)
gen l_b50noagg = ln(x50noagg)

gen dl_b50noagg = l_b50noagg - logb5575
gen dl_a50noagg = l_a50noagg - loga5575

reg l_a50noagg l_area 
predict l_a50noagg_narea, residual 
reg l_b50noagg l_area 
predict l_b50noagg_narea, residual 

reg loga5575 l_area 
predict loga5575_narea, residual 
reg logb5575 l_area 
predict logb5575_narea, residual 

#d;
twoway(scatter dl_b50noagg cbddist , sort color(black) msize(small))(lpoly logb5575_narea cbddist, sort lcolor(black) lwidth(medthick) lpattern(dash))(lfitci dl_b50noagg cbddist , sort acolor(gray%20) lcolor(black) lpattern(solid) lwidth(medthick)), 
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	legend(order(1 "Amentiies in 1950 (when no agglomeration) relative to average fundamental amenities in 1955-75"  2 "Polynomial fit of average fundamental amenities in 1955-75 (area size adjusted)" 4 "Fitted line of amenities in 1950 (when no agglomeration): slope = -0.06***") row(3) size(4))
	xtitle("Distance to CBD (kilometers)", size(5)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Log fundamental amenities in 1950 (without agglomeration forces)" "compared to the average between 1955-75 (dashed line)", size(4)) ylabel(-0.6 -0.4 -0.2 0 0.2 0.4, angle(horizontal) labsize(large)) yscale(r(-0.6 0.4));
	graph export "../../output/figure/figure_amenities_1950noagg_cbddist.pdf", as(pdf)  replace ;
#d cr

#d;
twoway(scatter dl_a50noagg cbddist , sort color(black) msize(small))(lpoly loga5575_narea cbddist, sort lcolor(black) lwidth(medthick) lpattern(dash))(lfitci dl_a50noagg cbddist, sort acolor(gray%20) lcolor(black) lpattern(solid) lwidth(medthick)), 
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	legend(order(1 "Productivity in 1950 (when no agglomeration) relative to average fundamental productivity in 1955-75" 2 "Polynomial fit of average fundamental productivity in 1955-75 (area size adjusted)" 4 "Fitted line of productivity in 1950 (when no agglomeration): slope = -0.07***") row(3) size(4))
	xtitle("Distance to CBD (kilometers)", size(5)) xlabel(, angle(horizontal) labsize(large))
	ytitle("Log fundamental productivity in 1950 (without agglomeration forces)" "compared to the average between 1955-75 (dashed line)", size(4)) ylabel(-0.6 -0.4 -0.2 0 0.2 0.4, angle(horizontal) labsize(large)) yscale(r(-0.6 0.4));
	graph export "../../output/figure/figure_productivity_1950noagg_cbddist.pdf", as(pdf)  replace ;
#d cr

* ***************************************************************
* FIGURE G.1 IN APPENDIX -- EVALUATION_NOISE_FIGURE *************
* ***************************************************************
import delimited "../../output/quant/output_evaluations_1945_50_noise.csv", clear

gen lpopdens45 = log(pop45/area)
gen lpopdens50 = log(pop50/area)

gen lpopdens50_high = log(pop50_high95/area)
gen lpopdens50_low = log(pop50_low95/area)

gen lempdens45 = log(emp45/area)
gen lempdens50 = log(emp50/area)

gen lempdens50_high = log(emp50_high95/area)
gen lempdens50_low = log(emp50_low95/area)

* Log population density -- Data vs Simulation ** Add weight by area size for polyfit 
#d;
twoway (lpoly lpopdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(scatter lpopdens50 cbddist  , sort color(gray%70) msymbol(Sh))
(scatter lpopdens50_low cbddist, sort color(black) msymbol(circle) msize(small)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (data)" 3 "1950 (lower bound when we have noise in 1945)" ) row(1) size(medium) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log population density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) yscale(range(2 12)); 
	graph export "../../output/figure/figure_density_pop_noise_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr


#d;
twoway (lpoly lempdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(scatter lempdens50 cbddist  , sort color(gray%70) msymbol(Sh))
(scatter lempdens50_low cbddist, sort color(black) msymbol(circle) msize(small)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (data)" 3 "1950 (lower bound when we have noise in 1945)" ) row(1) size(medium) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log employment density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) yscale(range(2 12)); 
	graph export "../../output/figure/figure_density_emp_noise_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr

* ************************************************************************
* FIGURE G.2 and G.3 IN APPENDIX -- EVALUATION_ROBUSTNESS_FIGURE *****************
* ************************************************************************
import delimited "../../output/quant/output_evaluations_1945_50_robustness.csv", clear 
gen lpopdens45 = log(pop45/area)
gen lpopdens50 = log(pop50/area)
gen lpopdens50_base = log(pop50_base/area)
gen lpopdens50_smallrho = log(pop50_smallrho/area)
gen lpopdens50_largetheta = log(pop50_largetheta/area)
gen lpopdens50_largegammah = log(pop50_largegammah/area)

gen lempdens45 = log(emp45/area) 
gen lempdens50 = log(emp50/area)
gen lempdens50_base = log(emp50_base/area)
gen lempdens50_smallrho = log(emp50_smallrho/area)
gen lempdens50_largetheta = log(emp50_largetheta/area)
gen lempdens50_largegammah = log(emp50_largegammah/area)

#d;
twoway (scatter lpopdens50_base lpopdens50, sort msymbol(Sh) color(gray) msize(small))  
(scatter lpopdens50_smallrho lpopdens50 , sort msymbol(Oh) color(brack) msize(small)) 
(scatter lpopdens50_largetheta lpopdens50 , sort msymbol(Th) color(brack) msize(small)) 
(scatter lpopdens50_largegammah lpopdens50 , sort msymbol(Dh) color(brack) msize(small)),
	legend(order(1 "1950 (base prediction in Section 5.5)" 2 "1950 (small rho)" 3 "1950 (large theta)" 4 "1950 (large floor space cost share)") row(2) size(medlarge))
	aspectratio(1) scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.5)
	xtitle("Log population density (per square km) in 1950 (Data)", size(5.0)) xlabel(4 5 6 7 8 9 10, angle(horizontal) labsize(medlarge)) xscale(r(4 10))
	ytitle("Log population density (per square km) in 1950 (Model prediction)", size(5.0)) ylabel(4 5 6 7 8 9 10, angle(horizontal) labsize(medlarge)) yscale(r(4 10)); 
	graph export "../../output/figure/figure_density_pop_1950_robustness.pdf", as(pdf)  replace ;
#d cr


#d;
twoway (scatter lempdens50_base lempdens50, sort msymbol(Sh) color(gray) msize(small))  
(scatter lempdens50_smallrho lempdens50 , sort msymbol(Oh) color(brack) msize(small)) 
(scatter lempdens50_largetheta lempdens50 , sort msymbol(Th) color(brack) msize(small)) 
(scatter lempdens50_largegammah lempdens50 , sort msymbol(Dh) color(brack) msize(small)),
	legend(order(1 "1950 (base prediction in Section 5.5)" 2 "1950 (small rho)" 3 "1950 (large theta)" 4 "1950 (large floor space cost share)") row(2) size(medlarge))
	aspectratio(1) scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.5)
	xtitle("Log employment density (per square km) in 1950 (Data)", size(5.0)) xlabel(4 5 6 7 8 9 10, angle(horizontal) labsize(medlarge)) xscale(r(4 10))
	ytitle("Log employment density (per square km) in 1950 (Model prediction)", size(5.0)) ylabel(4 5 6 7 8 9 10, angle(horizontal) labsize(medlarge)) yscale(r(4 10)); 
	graph export "../../output/figure/figure_density_emp_1950_robustness.pdf", as(pdf)  replace ;
#d cr

*****************************************************************************
* FIGURE 6 IN PAPER -- NOAGG_FIGURE *****************************************
*****************************************************************************
import delimited "../../output/quant/output_no_agglomeration.csv", clear 
save temp, replace 

gen lpopdens45 = log(pop45/area)
gen lpopdens50 = log(pop50/area)
gen lempdens45 = log(emp45/area)
gen lempdens50 = log(emp50/area)

gen lpopdens50_noagg = log(pop50_noagg/area)
gen lempdens50_noagg = log(emp50_noagg/area)

* Log population density in 1950 -- Data vs Simulation ** Add weight by area size for polyfit 
#d;
twoway (lpoly lpopdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lpopdens50_noagg cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lpopdens50 cbddist [aw=area], sort color(gray) lwidth(thick) lpattern(dash))  
(scatter lpopdens50_noagg cbddist, sort color(black) ) 
(scatter lpopdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log population density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) yscale(range(2 12)); 
	graph export "../../output/figure/figure_density_pop_noagg_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr

* Log Employment density in 1950 -- Data vs Simulation ** Add weight by area size for polyfit
#d;
twoway (lpoly lempdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lempdens50_noagg cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lempdens50 cbddist [aw=area], sort color(gray) lwidth(thick) lpattern(dash))
(scatter lempdens50_noagg cbddist , sort color(black) ) 
(scatter lempdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log employment density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) ; 
	graph export "../../output/figure/figure_density_emp_noagg_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr

******************************************************************************
* FIGURE 7 IN PAPER -- COUNTERFACTUALPATH_FIGURE *****************************
******************************************************************************
import delimited "../../output/quant/output_counterfactualpath.csv",clear 
save temp, replace 

gen lpopdens45 = log(pop45/area)
gen lpopdens50 = log(pop50/area)
gen lempdens45 = log(emp45/area)
gen lempdens50 = log(emp50/area)

gen lpopdens50_alt = log(pop50_alt/area)
gen lempdens50_alt = log(emp50_alt/area)
gen cbd_name =  "Counterfactual CBD" if lpopdens50_alt > 8 

* Log population density in 1950 -- Data vs Simulation ** Add weight by area size for polyfit 
#d;
twoway (lpoly lpopdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lpopdens50_alt cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lpopdens50 cbddist [aw=area], sort color(gray) lwidth(thick) lpattern(dash))  
(scatter lpopdens50_alt cbddist, sort color(black) mlabel(cbd_name) mlabsize(5) mlabcolor(blue) mlabpos(12)) 
(scatter lpopdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log population density (per square km)", size(7.0)) ylabel(0 2 4 6 8 10 12, angle(horizontal) labsize(huge)) yscale(range(2 12)); 
	graph export "../../output/figure/figure_density_pop_rationalcf_scatter_with_areasizeweights.pdf", as(pdf)  replace ;
#d cr


* Log Employment density in 1950 -- Data vs Simulation ** Add weight by area size for polyfit
#d;
twoway (lpoly lempdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lempdens50_alt cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lempdens50 cbddist [aw=area], sort color(gray) lwidth(thick) lpattern(dash))
(scatter lempdens50_alt cbddist , sort color(black) mlabel(cbd_name) mlabsize(5) mlabcolor(blue) mlabpos(12)) 
(scatter lempdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log employment density (per square km)", size(7.0)) ylabel(0 2 4 6 8 10 12, angle(horizontal) labsize(huge)) ; 
	graph export "../../output/figure/figure_density_emp_rationalcf_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr

*************************************************************************
* FIGURE G.7 IN APPENDIX -- FIGURE_LANDOWNER ****************************
*************************************************************************
import delimited "../../output/quant/output_no_agglomeration_landowner.csv", clear 
keep geocode area cbddist pop45 pop50 pop50_noagg emp45 emp50 emp50_noagg q50 q50_noagg
save temp, replace 

gen lpopdens45 = log(pop45/area)
gen lpopdens50 = log(pop50/area)
gen lempdens45 = log(emp45/area)
gen lempdens50 = log(emp50/area)

gen lpopdens50_noagg = log(pop50_noagg/area)
gen lempdens50_noagg = log(emp50_noagg/area)

* Log population density -- Data vs Simulation ** Add weight by area size for polyfit 
#d;
twoway (lpoly lpopdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lpopdens50_noagg cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lpopdens50 cbddist [aw=area], sort color(gray) lwidth(thick) lpattern(dash))  
(scatter lpopdens50_noagg cbddist, sort color(black) ) 
(scatter lpopdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.55)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log population density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) yscale(range(2 12)); 
	graph export "../../output/figure/figure_density_pop_noagg_landowners_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr


* Log Employment density -- Data vs Simulation ** Add weight by area size for polyfit
#d;
twoway (lpoly lempdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lempdens50_noagg cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lempdens50 cbddist [aw=area], sort color(gray) lwidth(thick) lpattern(dash))
(scatter lempdens50_noagg cbddist , sort color(black) ) 
(scatter lempdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.55)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log employment density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) ; 
	graph export "../../output/figure/figure_density_emp_noagg_landowners_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr


*********************************************************************************
* FIGURE G.8 IN APPENDIX -- EVALUATION_MYOPIC_FIGURE ****************************
*********************************************************************************
import delimited "../../output/quant/output_evaluations_1945_50_myopic.csv", clear 

gen lpopdens45 = log(pop45/area)
gen lpopdens50 = log(pop50/area)
gen lpopdens50_p = log(pop50_memory/area)

gen lempdens45 = log(emp45/area) 
gen lempdens50 = log(emp50/area)
gen lempdens50_p = log(emp50_memory/area)

* Log population density -- Data vs Simulation ** Add weight by area size for polyfit 
#d;
twoway (lpoly lpopdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lpopdens50_p cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lpopdens50 cbddist [aw=area], sort color(gray)lwidth(thick) lpattern(dash))
(scatter lpopdens50_p cbddist, sort color(black)) 
(scatter lpopdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log population density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) yscale(range(2 12)); 
	graph export "../../output/figure/figure_density_pop_myopic_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr

* Log Employment density -- Data vs Simulation ** Add weight by area size for polyfit
#d;
twoway (lpoly lempdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lempdens50_p cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lempdens50 cbddist [aw=area], sort color(gray) lwidth(thick) lpattern(dash))
(scatter lempdens50_p cbddist , sort color(black) ) 
(scatter lempdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log employment density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) ; 
	graph export "../../output/figure/figure_density_emp_myopic_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr

*********************************************************************************
* FIGURE G.9 IN APPENDIX -- EVALUATION_MEMORY_FIGURE ****************************
********************************************************************************* 
import delimited "../../output/quant/output_evaluations_1945_50_memory.csv", clear 

gen lpopdens45 = log(pop45/area)
gen lpopdens50 = log(pop50/area)
gen lpopdens50_p = log(pop50_memory/area)

gen lempdens45 = log(emp45/area) 
gen lempdens50 = log(emp50/area)
gen lempdens50_p = log(emp50_memory/area)

* Log population density -- Data vs Simulation ** Add weight by area size for polyfit 
#d;
twoway (lpoly lpopdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lpopdens50_p cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lpopdens50 cbddist [aw=area], sort color(gray)lwidth(thick) lpattern(dash))
(scatter lpopdens50_p cbddist, sort color(black)) 
(scatter lpopdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log population density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) yscale(range(2 12)); 
	graph export "../../output/figure/figure_density_pop_memory_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr


* Log Employment density -- Data vs Simulation ** Add weight by area size for polyfit
#d;
twoway (lpoly lempdens45 cbddist [aw=area], sort color(cranberry) lwidth(thick) lpattern(shortdash))  
(lpoly lempdens50_p cbddist [aw=area], sort color(black) lwidth(thick)) 
(lpoly lempdens50 cbddist [aw=area], sort color(gray) lwidth(thick) lpattern(dash))
(scatter lempdens50_p cbddist , sort color(black) ) 
(scatter lempdens50 cbddist  , sort color(gray%70) msymbol(Sh)),
	legend(order(1 "1945 (after bombing)" 2 "1950 (prediction)" 3 "1950 (data)" 4 "1950 (prediction)" 5 "1950 (data)") row(2) size(vlarge) holes(4))
	scheme(s1color) plotregion(lwidth(none) ilwidth(none)) scale(0.6)
	xtitle("Distance to CBD (km)", size(7.0)) xlabel(, angle(horizontal) labsize(huge))
	ytitle("Log employment density (per square km)", size(7.0)) ylabel(2 4 6 8 10 12, angle(horizontal) labsize(huge)) ; 
	graph export "../../output/figure/figure_density_emp_memory_scatter_with_areasizeweight.pdf", as(pdf)  replace ;
#d cr



